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Abstract 

An optimal method to constrain the non-linearity parameter <7nl °f the local- 
type non-Gaussianity from CMB data is proposed. The optimal estimator of <?nl 
is computed in the real space, which significantly reduces the computational cost 
and enables us to apply the method to the actual data. This method provides a 
much tighter constraint on g^L than any other suboptimal methods. Applying our 
method to the WMAP 9-year data, for example, we obtain the constraint $nl — 
(-3.3 ±2.2) x 10 5 , which is a few times tighter than previous ones. We also make a 
forecast for PLANCK data by using the Fisher matrix analysis. 



1 Introduction 



One of the most fundamental questions in cosmology is what is the origin of the primordial 
fluctuations which seed the large scale structure in the observed Universe as well as the 
anisotropies in the cosmic microwave background (CMB). Various cosmological observa- 
tions consistently show that primordial fluctuations are adiabatic, nearly scale-invariant 
and Gaussian, which is consistent with a prediction of simple single-field slow-roll inflation 
models. On the other hand, there are a variety of models in which probability distribution 
of primordial perturbations can significantly derivate from Gaussian ones. 

Among an infinite number of possibilities for deviation from Gaussian distributions, 
we in this paper focus on the local- type non-Gaussianity [1], in which the non-Gaussian 
curvature perturbation is given as a function of its Gaussian part $g(^D only at the 
same point, i.e. 

= $ G (x) + /nl [<& g(^) 2 - (* g0?) 2 >] + 4nl$g0?) 3 + • ■ ■ , (1) 

where /nl and <?nl are called non-linearity parameters. This type of non-Gaussianity is of 
particular interest. In standard single-field slow-roll inflation models, amount of this type 
of non-Gaussianity is too small to be observed at least in the near future. On the other 
hand, a range of theoretical models based on the inflationary Universe, in which multiple 
degrees of freedom during inflation contribute to primordial perturbations, can generate 
a large local- type non-Gaussianity [2]. Among them, typical examples are the curvaton 
scenario [3, 4, 5] and the modulated reheating scenario [6, 7]. Therefore the local- type 
non-Gaussianity is a unique probe for these models which may manifest only at the very 
early Universe and very high energy scales. 

So far many attempts have been made to detect /nl- Current constraints come from 
largely two types of observations. One is the bispectrum of the temperature anisotropy in 
CMB. Current data from the WMAP 9-year observation gives a constraint —3 < /nl < 
77 at 95% confidence level (C.L.) [8]. #1 Another probe is the scale-dependent bias in 
correlation functions of massive objects. At present the best constraint from the scale- 
dependent bias gives —37 < /nl < 25 at 95% C.L. [13], where angular correlation functions 
of galaxies and quasars are used. 

While /nl, which parameterizes the leading-order non-Gaussian term in Eq.(l), have al- 
ready begun to be constrained by various data, there are higher-order terms, which remain 
to be explored more deeply. Regarding <7nl, which is the coefficient of the next-to- leading 
order term in Eq. (1), so far only a few groups [17, 18] have presented CMB constraints on 
it. #2 In Ref. [17], the authors present a constraint -7.4 x 10 5 < £f NL < 8.2 x 10 5 (95% C.L.) 

While primordial perturbations are consistent with adiabatic ones and non-Gaussianity in curvature 
perturbations is discussed in the most literature, some theoretical models predict isocurvature perturba- 
tions which can have a local- type non-Gaussianity, e.g., Refs. [9, 10]. Current constraints on isocurvature 
non-Gaussianitics are presented in Refs. [11. 12]. 

* 2 Snl can be also constrained from the scale-dependent biases. For current constraints, we refer to 
Ref. [13]. 
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from the WMAP 5-year data based on suboptimal estimators. On the other hand, in Ref. 
[18], an optimal estimator of #nl is proposed. However, the estimator is computationally 
demanding and several approximations are required in application to actual data, which 
makes constraints less tight. The resultant constraint is g^L = (1.6 ± 7.0) x 10 5 at 1 a. 
Thus, so far there has been no optimal constraints on #nl from CMB which take into ac- 
count the inhomogeneity in noise levels of surveys and sky cuts for removal of foregrounds. 
Meanwhile there are models with large $nl which may be observationally detected (See, 
e.g., Refs. [14, 15]). * 3 Therefore it is important to improve constraints on c/nl and enable 
an optimal estimation of it from CMB observations. 

In this paper, we discuss a method to derive an optimal constraint on g^ from CMB 
observations. In particular, we show that an optimal estimator of ^nl can be computed 
much more efficiently in the real space than in the harmonic space. Organization of this 
paper is as follows. In the next section, we discuss a CMB trispectrum generated from 
<7nl- in Section 3, we present a representation of the optimal estimator for g-^L i n the 
real space, which we compute in this paper. In Section 4, after describing details of our 
analysis, we present constraints on ^nl from the WMAP 9-year data. We also compare 
our results with a forecast based on the Fisher matrix analysis in Section 5. The final 
section is devoted to summary. 

Throughout this paper, we assume a concordance flat power-law ACDM model, and 
the cosmological parameters are fixed to the mean values from the WMAP 7-year data 
alone [19], 

(fi 6 , fi c , H , r, n B , As) = (0.0448, 0.220, 71, 0.088, 0.963, 2.43 x 10" 9 ). (2) 

Here, Qb and Q c are respectively the density parameters for baryon and CDM, Ho is the 
Hubble constant in units of km/sec/Mpc, r is the optical depth of reionization, and n s 
and A s are respectively the spectral index and amplitude of power spectrum of curvature 
perturbations at a reference scale K = 0.002Mpc~\ i.e. P$(k) = iS-^rK^)™ 3 " 1 - 

2 CMB trispectrum from nonzero $nl 

First let us consider correlation functions of primordial curvature perturbations $ in the 
local-type non-Gaussianity. A non-trivial effect of $nl arises in the connected part of the 
four-point correlation function of $(x), or its Fourier dual, the trispectrum. If / NL = 0, 
the connected trispectrum should be given by 

mh)^(k 2 )^(k 3 )^(h)) conn = 6g NL [P^h)P^k 2 )P^k 3 ) + (3 perms)] (2tt) 3 <5^ (^1234), 

where ki x ... in = k^ + ■ ■ ■ + k in . For simplicity, we in this paper are to constrain <7nl, 
assuming / NL to be zero. 

# 3 We refer to Ref. [16] as a review of theoretical models which predict large local-type non-Gaussianities. 
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Neglecting the secondary non-Gaussianities arising from the second- or higher-order 
cosmological perturbation theory, the harmonic coefficients of the CMB temperature 
anisotropy from primordial perturbations $ can be given as 

rt,(fc)*(fc)^(fc), (4) 

where gi(k) is the temperature transfer function. The angular power spectrum Ci of the 
temperature anisotropy, which is defined by (ai m ai> m >*) = Ci5u>5 mm ', can be given as 

Q = - f dkk 2 gi (k) gi (k)P(k). (5) 



7T _ 

We define a reduced CMB trispectrum tz^^,* 4 

mi ^{21713^37713 ^Iattia) conn — ^IxIiIzIa ^mi 771277137714' 

(6) 

where G m [ 2 ^ m3 m 4 is defined by 

= J dfiY hm M)Yi, m MY hm MYh m M)- (7) 
Given the trispectrum of of Eq. (3), Eq. (6) leads to 

6 9Nt / dr ^KWAMfl,W&M + (3pe rm s)), (8) 
where ai(r) and 0i(r) are 

a i( r ) = \ J ' dkk 2 gi (k) 3l (kr), (9) 
Pi(r) = ^J dkk 2 gi (k) 3l (kr)P(k). (10) 
For later convenience, we introduce a normalized trispectrum t; 1 z 2 z 3 z 4 by 



d gN L 

6 / drr 2 (a h {r)(3 l2 (r)(3 h (r)(3 h (r) + (3 perms)) , (IT 



from which Eq. (8) can be recast into 



# 4 Note that our definition of a reduced trispectrum, fijj 2 ; 3 i 4 , is different from the one in Ref. [20], 



3 Fast optimal estimator of $nl 

According to Refs. [21, 18], an optimal estimator of <?nl can be given by 

K ■ 

9NL - r , 

\-fiprim/9 NL =l 

where an angle bracket with subscript #nl = 1 indicates an ensemble average over non- 
Gaussian simulations with unit #nl- Here K wlia is a quartic statistics, which is given 
by 



-Kprim — ^ ] ^ ] Mxl-ihljd m^raami [ a lrm 1 a l 2 m2 a hm 3 a l4m 4 

h—l4 mi —1714 



24 

d 

^^iimi Q'livri'z (^£31713^/41714)0 "I - 3 0^21712)0 (^£37713 ^47714 )o] , (14) 



where the bracket with subscript 0, (-)o, indicates an ensemble average over Gaussian 
simulations. Here a; m is a harmonic coefficient obtained from observed (or simulated) 
data maps weighted by inverse- variance, 

aim = y^Chji'mMm', (15) 
I'm' 

where Ci m y m < = C lmjVw/ + N lm j,t m > is the total covariance with Ci m) i> m , and N lmiVm , being 
denoted as those of signal and noise, respectively. 

As shown in Refs. [21, 18], Eq. (14) is proportional to the connected part of trispectrum 
of Eq. (3), so that the estimator Eq. (13) is unbiased. We can also see that 



(Ap r ; m Ap r ; m ^Q 
(-Kprim)g NI 



(<?NL/0 - ^2 

im 

1 



(-Kprim)c 



(16) 



'SNL=1 



where we used the relation (if pr im-Kprim)o = (-Kprim)p NL =i- See Appendix A for the proof. 
Then we finally obtain 



(2/3/4 

im2m 3 m i 

l' r ..l^m{-m' 4 

V V V V 1 _1 



x C, ,1 1 C, ,1 1 C, ,1 1 C, ,/ , ti' ji 11 v Q 1 / 2 3 ; 4 / / . (17) 

l\m\,l\m\ (21772, t^m^ '3 m 3>'3 m 3 t4m4,t 4 m^ l i l 2 l 3 b 4 m^rn^m^m^ \ > 

The right hand side is the inverse of the Fisher matrix of Eq. (21) with inhomogeneous 
noise and sky cuts being taken into account. Thus Eq. (17) shows that Eq. (13) is a 
minimum- variance estimator in the limit of weak non-Gaussianity. 
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Computation of A prim and (ifp r im) flNI =i should be implemented in the real space rather 
than the harmonic space. Combined with Eqs. (7) and (11), Eq. (14) can be rewritten in 
the real space as 

A prim = j d 3 r [A(r)B(r) 3 - 3A(r)B(r)(B(r) 2 ) 

-3B(f) 2 (A(f)B(f)) + 3(A(f)B(r)) (B(f) 2 ) ] , (18) 
where A(r) and B(r) are defined as 

A(r) = J2 a i( r ) Y im(r)a lm , (19) 

Im 

B(7 f ) = J %2Pi(r)Y lm (r)a lm . (20) 

im 

The form of Eq. (18) is computationally as demanding as the fast estimator of /nl of 
Refs. [22, 23, 24]. 



4 Constraint from the WMAP 9-year data 

In order to determine the normalization (ifprim) ffNI =i, we need to simulate non-Gaussian 
CMB maps. We adopt the method for local-type non-Gaussianity developed in Ref. [25], 
which uses the Gauss quadrature method with optimized nodes and weights for line of 
sight integral. Since the method is performed in the real space rather than the harmonic 
one, it is straightforwardly extended to the cubic model, while it is originally developed 
for models with up to the quadratic term in Eq. (1). In our analysis, we ask the method 
to be accurate with mean square of the error less than 0.01 at each multipole (Im). For 
^max = 1024, we found that this level of accuracy requires 42 quadrature nodes. 

Computation of inverse- variance weighted map di m is performed based on the method 
of Ref. [26], which uses a conjugate gradient method with multi-grid preconditioning. 
This method also allows to marginalize over amplitudes of components when their spatial 
template maps are provided. 

The transfer function of CMB is computed using the CAMB code [27]. We combine 
the foreground- cleaned maps of V and W bands of the WMAP 9- year data [8].* 5 with a 
resolution N si( ^ c = 512 in the HEALPix pixelization scheme [28]. * 6 We adopt the KQ75y9 
mask [8] which cuts 31.2% of the sky. We also set the maximum multipole Z max to 1024 
in our analysis. We marginalize the amplitudes of the monopole I = and dipoles I = 
1 as default and also optionally marginalize the amplitudes of the Galactic foreground 
components at large angular scales using the synchrotron, free-free and dust emission 
templates from Ref. [8]. 

# 5 http: / /lambda. gsfc.nasa.gov 
# 6 http://hcalpix.jpl. nasa.gov 
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Now we are to present our constraints on g^L- Without template marginalization 
of Galactic foregrounds, we obtain — (—3.9 ± 2.2) x 10 5 at 1 a. With template 
marginalization, this changes to #nl = (— 3.3±2.2) x 10 5 . Comparing these two constraints, 
we see that effects of residual Galactic foregrounds are not significant. As we shall show 
in the next section, sizes of the errors are almost the same as the expectation from the 
Fisher matrix analysis, which would in part support the validity of our analysis. Having 
all these results, we conclude that current CMB data is consistent with Gaussianity. 

As has been emphasized, our constraints are the first optimal ones on g^L from CMB. 
Thanks to the optimality of our method, our constraints are about a few times tighter than 
ones in previous studies [17, 18], which also use data of CMB temperature anisotropy. On 
the other hand, g^L can be also constrained from the scale-dependent bias in correlation 
functions of massive objects. Our constraints are almost comparable to one in Ref. [13], 
which uses data of galaxies and quasars. However, there would be one advantage in use 
of CMB data. While there is a significant degeneracy between / NL and #nl in constraints 
from the scale-dependent bias as seen in e.g. Ref. [13], there should be no degeneracy 
between these two parameters from CMB data at the leading-order both in /nl and <?nl- 
This is because in the limit of Gaussianity, a covariance between a bispectrum and a 
trispectrum of CMB anisotropies should vanish. 



5 Fisher matrix analysis 

We here present a method to evaluate an expected constraint on ^nl based on the Fisher 
matrix analysis, and apply it to the WMAP and forthcoming PLANCK surveys. 

Analogously to the case of bispectrum, the Fisher matrix for <?nl should be approxi- 
mately given as 

f , (f, , , , G'lkkU ^2 

rp _ J sky \ ^ \ A \ L hhhUz?rri\m2rnzrnA) rni \ 

where C\ = C\ + iVj is the total power spectrum and / S k y is a sky coverage. An error of 
#nl should be given by 1/ VF. 

Eq. (21) is computationally expensive as we at least have to carry out quadruple 
summation over Vs. However, using the fact that the trispectrum tj 1 j 2 j 3 i 4 of Eq. (11) is in 
a separable form, Eq. (21) can be brought into a computationally cheaper expression as 
follows: 



F = 48tt 2 J dfi J drr 2 J dr' r' 2 [C aa (r, r' , fi)C^(r, r', y) (22) 
+3C Q/3 (r, r', y)C pa {r, r', //)] C pf) {r, r', fif , 
where C aa ' (r, r' , /i), with a and a' being either a or /3, is given by 

c w (r!r>) = E (2i±iMKM P , (/l) . (23) 
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Here Pi(fj) is the Legendre polynomial. Contrast to the quadruple sum over Vs of Eq. (21), 
Eq. (23) has only a triple integrals, computation of which is thus modest. Derivation of 
Eq. (23) is presented in Appendix B. 

We compute C\ using the CAMB code [27]. On the other hand, we approximate N\ by 
the Knox's formula [29], 



N l = ^FWHM^T exp 



1(1 + 1)%™ 

V 7 8 In 2 



(24) 



where #fwhm is the full width at half maximum of the Gaussian beam, and a? is the root 
mean square of the instrumental noise par pixel. The total Ni of a multi-frequency survey 
can be given by a quadrature sum of iVj of each frequency band. 

First, we estimate the expected error on (?nl from the WMAP survey. Survey param- 
eters we adopted for the WMAP 9-year V and W bands for #fwhm and <jt are listed in 
Table 1. Furthermore, we assume that / s k y = 0.69, which is consistent with the mask we 
adopt in Section 4. From the above setup, we obtain A^nl = 2.1 x 10 5 . The size of error 
is almost the same as one we obtained from actual WMAP 9-year data in Section 4, which 
supports the validity of our analysis. 



band 


V 


W 


#fwhm [arcmin] 


21.0 


13.2 


Ot [hiK] 


21 


31 



Table 1: Survey parameters for 9-year observation of the WMAP V and W bands [30]. 



Next, we also forecast a constraint from a future survey. With the survey parameters 
for PLANCK listed in Table 2 and / sky = 0.69, we find that PLANCK will constrain # NL 
with error A^nl = 6.7 x 10 4 , which is a few times tighter than the current ones. 



band 


#fwhm [arcmin] 


33.0 


24.0 


14.0 


10.0 


7.1 


5.0 


5.0 


a T [mK] 


5.5 


7.4 


12.8 


6.8 


6.0 


13.1 


40.1 



Table 2: Same as in Table 1 but for the PLANCK 14-month observation [31]. 



6 Conclusion 

We present a method to derive an optimal constraint on the non-linearity parameter <?nl 
of the local-type non-Gaussianity from CMB data. Computational cost of our method is 
almost the same as the fast /nl estimator of the local-type non-Gaussianity. Applying 
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our method to the WMAP 9-year data, we obtain (?nl = (—3.3 ± 0.8) x 10 5 with tem- 
plate marginalization of Galactic foregrounds. The size of the error is consistent with an 
expectation based on the Fisher matrix analysis. 

Our constraints are a few times tighter than ones in the previous studies [17, 18] from 
CMB data. While the improvement is not very dramatic, our constraints are however the 
first optimal ones on #nl from CMB. We emphasize that by transferring from the harmonic 
space to the real space, we significantly reduced the computational cost in computing the 
optimal estimator. We expect our method can be applied to high resolution data of the 
PLANCK survey, and (?nl will be constrained tighter. 

In this paper we focused on the cubic term in Eq. (1). However, higher-order terms 
in the equation can also be nonzero. We insist that our method can be straightforwardly 
extended for any of the higher-order Gaussian terms of the local-type non-Gaussianity as 
far as we assume that all the other non-Gaussian terms should vanish. 
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A Equivalence of (K wim K wim ) and (i^ prim )^ NL =i 

In this appendix, we present a proof for (K pT i m K pT im)o = (-Kprim) 9NL =i- First, we divide 
the CMB anisotropy a\ m into the Gaussian a G> i m and non-Gaussian a^cim parts: 

a G ,im = 4:ir(-i) 1 

ClNG,lm = 47r(-i) Z 

which leads to a>i m = a>G,i m + ^nl^ng,;™- Since filtered map di m is a linear function of 
a>im, o,im can also be divided into the Gaussian and non-Gaussian parts in the same way, 



d 3 k 

(27p 



gi (k)<f> G (k)Y lm (k), 



(25) 
(26) 
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First we compute (iC pr i m ) 9NL =i. Up to the leading order in onc,;™, Eq. (14) leads to 



24 

h-'-li mi-- -1714 



■12{aNG,lim i a>G,l2m2){ a G,l 3 m 3 aG,l4m 4 )} 
777 ^ ^ ^ j thl2l3U@m1m2m.3m4 (27) 



24 

Zl-k "11—7714 

Z^-Z 4 mi -mi 



/' /' /' /' 

'1 '2*3*4 



where we used the relation 

(az imi a/ 2m2 a; 3m3 a/ 4m4 ) coim = g NL [{(aNG,Zimi5G,Z2m2^G,Z3m 3 5G,Um4) + (3 perms)} (28) 

-{(aNG,/im.iaG,«2m2)(aG,«3m3aG,/ 4 m 4 ) + (H perms)}] , 

and 

On the other hand, (if pr im-Kprim)o is a ls° computed from Eq. (14), which leads to 

(^prim-frprim)o = 7j~7^ ^k^mi^msm (^0) 

x ( SG,Z 1 miSG,J2OT.2^G,Z3m3 5G,i477i4 — 6aG,Z 1 miOG,i 2 m2 (aG,; 3 m 3 aG,Z4 m4 ) 
+3(a G ,Zimi ^G ; i277l2 )(aG,Z 3 m 3 aG,Z4m 4 ) 

By using the Wick theorem, after lengthy but simple calculation, we obtain 

(-Kprimi^prim)o — ^ ^ ^ ^ ^hhhU^ mim\mzm4 

(31) 

ti-i 4 mi-m 4 
J/...J/ mi-m 4 



^G^imi a G,Z 2 m 2 ^G,Z 3 m 3 ^G,Z^m 4 



' ;' ;' /' 

l'2 l 3*. 



^ ZlTTll ,ii77lj ^Z2'"2i'2 m 2^'3 m 3»'3 m 3^ 4m4 ''4 m 4^1^'3' 4 ^"l / 1 

Comparison of Eqs. (27) and (31) shows (K prim K prim } = (K pvim ) gNh=1 . 



4 

' 2 m 3 m 4 - 
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B Derivation of Fisher matrix 

In this appendix, we derive Eq. (23). We first start from computation of the summation 
over m's in Eq. (21). Using Eq. (7), we obtain 



4 

J2(&^ m J 2 = JdnJ dn'X\^Y hm Xn)Y* mi {n')] (32) 



{m} i=l rn,i 

_ (2/ 1 + l)(2/ 2 + l)(2/ 3 + l)(2/ 4 + l 
32tt 2 

where in the second equality we used the relation 

91 4- 1 

Y. Y ^) Y Lin') = ■ n'). (33) 

m 

Then adopting the definition of ti^i^ of Eq. (11), Eq. (21) can be rewritten as 

'(2/ + l)af ) (r)a; W (r') 



F= 12it 2 J drr 2 J dr' r' 2 J d^\\ 

a, a' i=l 



AtiCi 

d.d' t=l . 



(34) 



where a\ (r) and a ; (r) should be either ai(r) and 0i(r). We also introduced a vector 
a = (a^\ a^ 2 \ a^ 3 \ a^), which symbolically represents a permutation of (a, j3, j3, j3). By 
taking summations over a and b of any possible permutations, then we finally obtain Eq. 
(23). 
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